home *** CD-ROM | disk | FTP | other *** search
/ MacWorld 1997 September / Macworld (1997-09).dmg / Serious Software / Cherwell Scientific Demos / pro Fit / pro Fit 5.0 demo (fpu).sea / pro Fit 5.0 demo (fpu) / External Modules / External modules sources / C / FFT.c next >
Text File  |  1996-04-28  |  11KB  |  329 lines
  1. /*
  2.  * This unit contains the function 'fft' which, together with 350 additional
  3.  * programs useful in scientific work, is found in the book 'Numerical Recipes:
  4.  * The Art of Scientific Computing', available from Cambridge University Press.
  5.  * It is used here by permission of the authors.
  6.  */
  7.  
  8.  
  9. /***************************************************************************************/
  10. /* FFT.c                                                                                  */
  11. /*                                                                                     */
  12. /* Version 26.9.94                                                                     */
  13. /***************************************************************************************/
  14.  
  15. #include "proFit_interface.h"
  16. #include <math.h>
  17.  
  18.  
  19.  
  20. static double sqr(double x) {return (x*x);}
  21.  
  22. static void fft (double* x, int num,char inverse)
  23.     {
  24.         long ii, jj, n, mmax, m, j, istep, i;
  25.         double wtemp, wr, wpr, wpi, wi, theta, tempr, tempi;
  26.     
  27.         n = num;
  28.         j = 1;
  29.         for (ii = 1;ii<= (n / 2);ii++) 
  30.             {
  31.                 i = 2 * ii - 1;
  32.                 if (j > i) 
  33.                     {
  34.                         tempr = x[j];
  35.                         tempi = x[j + 1];
  36.                         x[j] = x[i];
  37.                         x[j + 1] = x[i + 1];
  38.                         x[i] = tempr;
  39.                         x[i + 1] = tempi;
  40.                     }
  41.                 m = n / 2;
  42.                 while ((m >= 2) && (j > m)) 
  43.                     {
  44.                         j = j - m;
  45.                         m = m / 2;
  46.                     }
  47.                 j = j + m;
  48.                 if (TestStop())         // test if proFit or a user tries to stop the program }
  49.                     goto end;
  50.             }        //for ii = 1 to (n div 2) do }
  51.  
  52.         mmax = 2;
  53.         while (n > mmax) 
  54.             {
  55.                 istep = 2 * mmax;
  56.                 theta = 6.283185307959 / mmax;
  57.                 if (inverse) 
  58.                     theta = -theta;
  59.                 wpr = -2 * sqr(sin(0.5 * theta));
  60.                 wpi = sin(theta);
  61.                 wr = 1;
  62.                 wi = 0;
  63.                 for (ii = 1;ii<= (mmax / 2);ii++) 
  64.                     {
  65.                         m = 2 * ii - 1;
  66.                         for (jj = 0;jj<= ((n - m) / istep);jj++) 
  67.                             {
  68.                                 i = m + jj * istep;
  69.                                 j = i + mmax;
  70.                                 tempr = wr * x[j] - wi * x[j + 1];
  71.                                 tempi = wr * x[j + 1] + wi * x[j];
  72.                                 x[j] = x[i] - tempr;
  73.                                 x[j + 1] = x[i + 1] - tempi;
  74.                                 x[i] = x[i] + tempr;
  75.                                 x[i + 1] = x[i + 1] + tempi;
  76.                             }        // for jj := 0 to ((n - m) div istep) do }
  77.                         wtemp = wr;
  78.                         wr = wr + wr * wpr - wi * wpi;
  79.                         wi = wi + wi * wpr + wtemp * wpi;
  80.                         if (TestStop())     //test if proFit or a user tries to stop the program }
  81.                             goto end;
  82.                     }    // for ii := 1 to (mmax div 2) do }
  83.                 mmax = istep;
  84.             }        // while n>mmax }
  85.             end:;
  86.     }        // fft }
  87.  
  88.  
  89. static    void shiftFFT (double* x, long np)
  90.     {
  91.         double xx;
  92.         long i;
  93.         for (i = 1;i<= np;i++) 
  94.             {
  95.                 xx = x[i];
  96.                 x[i] = x[i + np];
  97.                 x[i + np] = xx;
  98.             }
  99.     }    // shiftFFT }
  100.  
  101. static void CleanUpRun(double* x)
  102. {    StopExecution();
  103.     if (x!=0)  {SysBeep(1);DisposePtr((Ptr) x);x=0;}
  104. }
  105.  
  106.  
  107.  
  108. /***************************************************************************************/
  109.  
  110. void SetUp (    short* const moduleKind,        /* set moduleKind to isFunction or isProgram */
  111.                 Str255 name,                    /* the name of the program or function (pascal string) */
  112.                 long* const requiredGlobals,    /* the number of bytes to be allocated in ExtModulesParamBlock.globals */
  113.                                                 /* set requiredGlobals to 0 if you don't use this feature */
  114.                 ExtModulesParamBlock* pb)        /* the complete parameter block passed by pro Fit to the */
  115.                                                 /* routines defined in this file. In most cases it can be ignored */
  116. /* SetUp is called once when the external module is linked to proFit */
  117. {
  118.     *moduleKind=isProgram;                        /* we define a program */
  119.     SetPascalStr(name,"\pComplex FFT",255);        /* with the name "Complex FFT" */
  120.     *requiredGlobals=0;                            /* we define no globals */
  121. }
  122.  
  123.  
  124. /***************************************************************************************/
  125.  
  126. void InitializeProg (ExtModulesParamBlock* pb)
  127.     /* Can be left emtpy if not needed. */
  128.     /* called when the external module is linked to proFit after SetUp was called */
  129.     /* can be used to inititialize global variables, etc. */
  130. {        pb->v[0] = 0;    //use it as boolean value to test for the first time this program runs }
  131.         pb->v[1] = 64;    //these are the default values for the dialog box asking for dataNr, ColumnNrs 
  132.         pb->v[2] = 1;
  133.         pb->v[3] = 4;
  134.         pb->v[4] = 0;
  135.         pb->v[5] = 0;
  136. }
  137.  
  138.  
  139.  
  140. /***************************************************************************************/
  141.  
  142. void Run(ExtModulesParamBlock* pb)
  143. /* pro Fit calls this function when the name of the program is chosen from the */
  144. /* Run Program submenu in the menu Calc */
  145. {            InputRec r;
  146.             double    ex1=pb->v[1],
  147.                     ex2=pb->v[2],
  148.                     ex3=pb->v[3],
  149.                     ex4=pb->v[4],
  150.                     ex5=pb->v[5],
  151.                     x0, xinc, y0, yinc, yoff;
  152.             long i, xCol, yCol;
  153.             double* x=0;        // pointer to the data array }
  154.             long nrData;        // the number of data (integer power of 2) }
  155.             char inverse;    // false if normal FFT should be done, true for reverse FFT }
  156.             short answer;
  157.         
  158.         if (pb->v[0] == 0)
  159.             {        // this is the first time this program runs }
  160.                 pb->v[0] = 1;
  161.             //\p instructs the compiler to produce a pascal string
  162.                 if (AskBox(&answer,"\pDo you want to get some information on this program in the Results window ?"))
  163.                     return;                // stop button pressed 
  164.                 else 
  165.                   {if (answer == 1)     //yes button pressed }
  166.                     {
  167.                         Writeln("\pFast Fourier Transform: converts between data.");
  168.                         Writeln("\pin the 'time domain' and data in the 'frequency domain'");
  169.                         Writeln("\p");
  170.                         Writeln("\pThe data in the time domain is arranged in");
  171.                         Writeln("\pthree columns:");
  172.                         Writeln("\p  1st column: the time for each data point");
  173.                         Writeln("\p  2nd column: the real value of each data point");
  174.                         Writeln("\p  2rd column: the imaginary value of each data point");
  175.                         Writeln("\pThe data in the time domain is also arranged in");
  176.                         Writeln("\pthree columns:");
  177.                         Writeln("\p  1st column: the frequency for each data point");
  178.                         Writeln("\p  2nd column: the real value of each data point");
  179.                         Writeln("\p  2rd column: the imaginary of value each data point");
  180.                         Writeln("\p");
  181.                         Writeln("\pWhen running the program, you must specify the");
  182.                         Writeln("\pnumber of data points (between 4 and 4096).");
  183.                         Writeln("\pThe number of data points should be a power of 2.");
  184.                         Writeln("\pIf it is not a power of 2, it is rounded to the");
  185.                         Writeln("\pnext power of 2.");
  186.                         Writeln("\pIf you run an FFT, set 'First input column'");
  187.                         Writeln("\pto the time column of the data points in the");
  188.                         Writeln("\ptime domain, 'First output column' to the");
  189.                         Writeln("\pfrequency column of the data points in the");
  190.                         Writeln("\pfrequency domain.");
  191.                         Writeln("\pIf you run an inverse FFT, set 'First input column'");
  192.                         Writeln("\pto the frequency column of the data points in the");
  193.                         Writeln("\pfrequency domain, 'First output column' to the");
  194.                         Writeln("\ptime column of the data points in the");
  195.                         Writeln("\ptime domain.");
  196.                         Writeln("\p'Offset of FFT' defines an additive constant for");
  197.                         Writeln("\pthe frequency values.");
  198.                         return;
  199.                     }
  200.                   }
  201.             }
  202.  
  203.         r[0].x =&ex1;
  204.         r[0].s = "\pNumber of data points";
  205.         r[1].x = &ex2;
  206.         r[1].s = "\p$CFirst input column";
  207.         r[2].x = &ex3;
  208.         r[2].s = "\p$CFirst output column";
  209.         r[3].x = &ex4;
  210.         r[3].s = "\pOffset of FFT";
  211.         r[4].x = &ex5;
  212.         r[4].s = "\p$XInverse FFT";
  213.         if (InputBox(5, &r)) return;    // if something is not ok we quit 
  214.         else
  215.             {    // if everything is ok we do something 
  216.  
  217.                 ex1 = fabs(ex1);
  218.                 if ((ex1 > 4096) || (ex1 < 4)) return;
  219.                 i = ex1;
  220.                 nrData = 4;
  221.                 while (2*nrData<=i) nrData *= 2;
  222.                 ex2 = fabs(ex2);
  223.                 if ((ex2 > 97) || (ex2 < 1))  return;
  224.                 xCol = (ex2);
  225.                 ex3 = fabs(ex3);
  226.                 if ((ex3 > 97) || (ex3 < 1))  return;
  227.                 yCol = ex3;
  228.                 yoff = ex4;
  229.                 inverse = (ex5 != 0);
  230.                 pb->v[1] = nrData;    //these are the default values for the dialog box asking for dataNr, ColumnNrs 
  231.                 pb->v[2] = ex2;
  232.                 pb->v[3] = ex3;
  233.                 pb->v[4] = ex4;
  234.                 pb->v[5] = inverse;
  235.  
  236.                 x = (double*) NewPtr(sizeof(double)*(2L*nrData+1));            // prepare memory for data 
  237.                 if (x == 0) return;
  238.  
  239.  
  240.                 if (TestData(1,xCol)) 
  241.                     {
  242.                         x0 = GetData(1,xCol);
  243.                         if (TestData(2,xCol))
  244.                             xinc =GetData(2,xCol) - x0;
  245.                     }
  246.  
  247.                 
  248.                 for (i = 1;i<= nrData;i++)            // read data from window 
  249.                     {
  250.                         if (TestData(i,xCol + 1))    //test if data value is ok }
  251.                             x[2 * i - 1] = GetData(i,xCol + 1);
  252.                         else
  253.                             x[2 * i - 1] = 0;
  254.                         if (TestData(i,xCol + 2))
  255.                             x[2 * i] = GetData(i,xCol + 2);
  256.                         else
  257.                             x[2 * i] = 0;
  258.                         if (TestStop()) {CleanUpRun(x);return;}
  259.                     }
  260.                 
  261.                 fft(x, 2 * nrData, inverse);    // fft calculations }
  262.                 yinc = 1 / (xinc*nrData) ;
  263.                 if (!inverse) 
  264.                     {
  265.                         shiftFFT(x, nrData);
  266.                         y0 = yoff - 1 / xinc / 2;
  267.                     }
  268.                 else
  269.                     {
  270.                         y0 = yoff;
  271.                         for (i = 1;i<= nrData / 2;i++) 
  272.                             {
  273.                                 x[4 * i - 3] = x[4 * i - 3] / nrData;
  274.                                 x[4 * i - 2] = x[4 * i - 2] / nrData;
  275.                                 x[4 * i - 1] = -x[4 * i - 1] / nrData;
  276.                                 x[4 * i] = -x[4 * i] / nrData;
  277.                             }
  278.                     }
  279.                 if (TestStop() ) {CleanUpRun(x);return;}
  280.                 SetColName(yCol,"\pfrequency");
  281.                 SetColName(yCol + 1,"\pFFT real");
  282.                 SetColName(yCol + 2, "\pFFT imag.");
  283.                 for (i = 1; i<= nrData;i++)            // output of data into window }
  284.                     {
  285.                         SetData(i,yCol, y0 + (i - 1.0) * yinc);
  286.                         SetData(i,yCol + 1, x[2 * i - 1]);
  287.                         SetData(i,yCol + 2, x[2 * i]);
  288.                         if (TestStop() )  {CleanUpRun(x);return;}
  289.                     }
  290.  
  291.                 if (x != 0) 
  292.                 {    DisposPtr((Ptr) x);
  293.                     x=0;
  294.                 }        // we don't need the memory any more }
  295.             }
  296. }
  297.  
  298. /***************************************************************************************/
  299.  
  300. void CleanUp (ExtModulesParamBlock* pb)
  301.     /* called when the function or program is removed from pro Fit's menus */
  302.     /* in most cases, this function can be empty */
  303. {
  304. }
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  
  313. /***************************************************************************************/
  314.                         /* for functions, not used here: */
  315. /***************************************************************************************/
  316.  
  317. void InitializeFunc (Boolean* const hasDerivatives, Str255 descr1stLine, Str255 descr2ndLine,        
  318.                     short* const numberOfParams, DefaultParamInfo* const a0, ExtModulesParamBlock* pb)
  319. {}
  320. void Func (    double x, ParamArray a,    double* const y, ExtModulesParamBlock* pb)        
  321. {}
  322. void Derivatives(double x, ParamArray a, ParamArray dyda, ExtModulesParamBlock* pb)
  323. {}
  324. short Check(short paramNo, DefaultParamInfo* const a0, ExtModulesParamBlock* pb)
  325. {return ok;}
  326. void First (ParamArray a, ExtModulesParamBlock* pb)
  327. {}
  328. void Last (ExtModulesParamBlock* pb)
  329. {}